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The interaction between the supermassive black hole at the centre of the 
Milky Way, Sagittarius A’, and its accretion disk occasionally produces 


high-energy flares seen in X-ray, infrared and radio. One proposed 
mechanism that produces flares is the formation of compact, bright 
regions that appear within the accretion disk and close to the event horizon. 
Understanding these flares provides a window into accretion processes. 
Although sophisticated simulations predict the formation of these flares, 
their structure has yet to be recovered by observations. Here we show a 
three-dimensional reconstruction of an emission flare recovered from 
Atacama Large Millimeter/Submillimeter Array light curves observed on 
11 April 2017. Our recovery shows compact, bright regions at a distance 

of roughly six times the event horizon. Moreover, it suggests a clockwise 
rotation in alow-inclination orbital plane, consistent with prior studies 
by GRAVITY and the Event Horizon Telescope. To recover this emission 


structure, we solve an ill-posed tomography problem by integrating a neural 
three-dimensional representation with a gravitational model for black holes. 
Although the recovery is subject to, and sometimes sensitive to, the model 
assumptions, under physically motivated choices, our results are stable and 


our approach is successful on simulated data. 


The compact region around the Galactic Centre supermassive black 
hole Sagittarius (Sgr) A’ isa unique environment where the magnetized 
turbulent flow of an accretion disk is subject to extreme gravitational 
physics. The dynamical evolution of this complex system occasionally 
leads to the production of energetic flares’ seen in X-ray’, infrared’ and 
radio*. The physical nature, structure, origin, formation and even- 
tual dissipation of flares are topics of active research*” * key to our 
understanding of accretion flows around black holes. One proposed 
explanation for Sgr A flares is the formation of compact bright regions 
caused by hot pockets of lower-density plasma within the accretion 
disk, which are rapidly energized (for example, through magnetic 
reconnection’). These ‘bubbles’, ‘hotspots’ or ‘flux tubes’ observed in 


numerical simulations (for example, ref. 10) are hypothesized to form 
in orbit close to the innermost stable circular orbit (ISCO) of Sgr A’. The 
association of flares with orbiting hotspots close to the event hori- 
zonis consistent with near-infrared detections made by the GRAVITY 
Collaboration” and radio observations of the Atacama Large 
Millimeter/Submillimeter Array (ALMA). 

The context for this work is set by the first images" of Sgr A’ 
revealed by the Event Horizon Telescope (EHT) collaboration. The 
images, reconstructed from very-long-baseline interferometry obser- 
vations from 6-7 April 2017, show a ring-like structure with a central 
brightness depression—a strong suggestion that the source is indeed 
asupermassive black hole”. Even in its quiescent state, imaged by EHT 
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Fig. 1 | A 3D recovery ofa Sgr A flare observed by ALMA on 11 April 2017. a, The 
validation-y’, a robust data-fitting metric (Methods), indicates a preference of 
low inclination angles, 8, < 18°, with a local minimum around @, = 12° (red curve). 
For each inclination, the 3D recovery is run with five random initializations, 
producing a spread that indicates recovery stability. The blue curve indicates 
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that the analysis is largely insensitive to the black-hole spin. b, A recovered 3D 
volume visualized from two view angles in intrinsic (flat space) coordinates 
(the event horizon illustrated for size comparison). The recovery shows two 
emission regions (blue arrows) at radii of 11-13M (approximately six times the 
Schwarzschild radius). 


on 6-7 April, Sgr A’ has shown considerable structural variability. On 
11 April 2017, Sgr A' was observed by ALMA directly after a high-energy 
flare seenin X-ray. The ALMA light curves exhibit an even higher degree 
of variability than 6-7 April*”, including distinct coherent patternsin 
the linear polarization” with variability on the scale of an orbit. The 
presence of synchrotron-radiating matter very close to the horizon of 
Sgr A' could potentially give rise to bright three-dimensional (3D) struc- 
tures that orbit and evolve within the accretion disk. In this work, we 
present a 3D recovery of emission in orbit around Sgr A’, reconstructed 
from ALMA light curves observed on 11 April 2017 (Fig. 1). 

To achieve this 3D reconstruction, we developed anew computa- 
tional approach that we call orbital polarimetric tomography. In con- 
trast to prior work by refs. 11,13, which employed a strongly constrained 
parametric hotspot model with only a handful of parameters to tune 
and interpret the observations, the goal of this work is to recover the 
complex 3D structure of flares as they orbit and evolve in the accretion 
disk around Sgr A’. 

Tackling this inverse problem necessitates a change from typical 
tomography, wherein 3D recovery is enabled by multiple viewpoints. 
Instead, the tomography setting we propose relies on observing a 
structure in orbit, travelling through curved space-time, from a fixed 
viewpoint. As it orbits the black hole, the emission structure is observed 
(projected) along different curved ray paths. These observations of the 
evolving structure over time effectively replace the observations from 
multiple viewpoints required in traditional tomography. Our approach 
builds upon prior work on dynamical imaging and 3D tomography in 
curved space-time, which showed promising results in simulated future 
EHT observations”. 

Similar to the computational images recovered by EHT“, our 
approach solves an underconstrained inverse problem to fit amodelto 
the data. Nevertheless, ALMA observations do not resolve event horizon 
scales (-10° lower resolution), which makes the tomography problem we 
propose particularly challenging. To put it differently, we seek to recover 
an evolving 3D structure from a single-pixel observation over time. To 
solve this challenging task, we integrate the emerging approach of neural 
3D representations”, which has an implicit regularization that favours 
smooth structures” with physics constraints (details in Methods). The 
robustness of the results thus relies on the validity of the constraints 
imposed by the gravitational and synchrotron emission models. 

We take advantage of the very high signal-to-noise and cadence 
of the ALMA dataset’, as well as the linear polarization information”. 
The choice to only fit the linear polarization (LP) light curves reflects 
the uncertainty associated with the unpolarized intensity of the back- 
ground accretion disk. Although the total intensity light curve is domi- 
nated by the accretion disk, such extended emission structures are 


partially depolarized in animage-average sense”. In contrast, compact 
bright sources, suchas a putative hotspot, are characterized by a large 
fractional LP and fast evolution on dynamical timescales’, hence 
allowing separation of the flare component from the background 
accretion. In Supplementary Information Section 2.2, we quantita- 
tively assess the effect of the background accretion disk on simulated 
reconstruction results. 


Results 
On 11 April 2017, ALMA observed Ser A at -230 GHzas part of a larger 
EHT campaign (Fig. 2, top). The radio observations directly followed a 
flare seen in the X-ray. The LP, measured by ALMA-only light curves*”* 
as a complex time series Q(t) + iU(t), appears to evolve ina structured, 
periodic manner suggesting a compact emission structure in orbit. 
The work of ref. 13 hypothesizes a simple bright spot (that is, idealized 
point-source” or spherical Gaussian”) at radius (r) ~ 11M (where Mis the 
black-hole mass; 2M is the Schwarzschild radius); however, a rigorous 
data fitting was not performed. Furthermore, the proposed parametric 
model is limited and does not explain all the data features. The orbital 
polarimetric tomography approach that we propose enables a rigorous 
data fitting and recovery of flexible 3D distributions of the emitting 
matter, relaxing the assumption of a fixed orbiting feature enforced 
by prior studies’”*”°. This opens a new window into understanding the 
spatial structure and location of flares relative to the event horizon. 
Our model, detailed in Methods, is able to fit the ALMA light curve 
data very accurately (Fig. 2, bottom). The optimization procedure 
simultaneously constrains the inclination angle of the observer and 
estimates a 3D distribution of the emitting matter associated with this 
flaring event, starting from 9:20 UT (~30 min after the peak of the X-ray 
flare”). Despite the fact that ALMA observations are unresolved (effec- 
tively a single pixel with time-dependent complex LP information) at 
the horizon scale, our analysis suggests some interesting insights: 


e Lowinclination angles (0, < 18°, Fig. 1a, red) are preferred by the 
validation-x* (Methods). Although the methodology is different, 
this result is broadly consistent with EHT findings from 6-7 April”, 
which favoured low inclination angles of -30° by comparing recov- 
ered images with general relativistic (GR) magnetohydrodynamic 
simulations. The fiducial model of ref. 13 corresponded to an 
inclination angle of ~22°. Lowinclination was also favoured inthe 
analysis of the GRAVITY infrared flares” >”. 

¢ The recovered 3D emission has two compact bright regions at 
r=11M and 13M (Fig. 1b). The location (radius and azimuthal 
position) of the bright region is consistent with the qualitative 
analysis of refs. 13,26. 
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Fig. 2| ALMA light curves and a model fit over a period of -100 min. Top: 

229 GHz light curves were observed on 11 April 2017 (MJD 57854) as part of the 
EHT Sgr A’ campaign. The red-shaded region corresponds to atime period of 
~100 min in which polarimetric (Q-U) loops are apparent, directly after an X-ray 
flare was observed (grey-shaded region). The rotation of the polarization angle at 
a period similar to a Keplerian orbital period suggests the signal is coming froma 


bright compact structure in orbit around Sgr A”. Bottom: a data fit is preformed 
onthe intrinsic LP curves (centred and derotated). The model light curves are 
produced through ray tracing the estimated 3D volume at a fiducial inclination 
angle of 8, = 12°. The resulting model light curves accurately describe the data, 
including the small looping feature highlighted by the blue arrow. 


Data fitting 

Before solving the tomography problem, we perform preprocessing 
according to the procedure outlined in ref. 13. In particular, we time 
average the data, subtract a constant (time-averaged) LP component 
interpreted as the ring-like accretion disk component observed by 
the EHT and derotate the electric vector polarization angle to account 
for Faraday rotation (details in Methods). Figure 2 illustrates the data 
before and after the preprocessing. 

To obtain a model prediction, a3D emission structure is adjusted 
so that when placed in orbit, the numerically ray-traced light curves 
align with the observations. To recover the vertical structure, our 
approach primarily leverages asymmetries in the polarimetric radiative 
transfer. In particular, the geometry of space-time and the magnetic 
field dictate the angle of linear polarization (Q-U). Moving an emission 
point changes the observed angle of linear polarization. Thus, errone- 
ously placing emission at time t = O and propagating it in time will rotate 
to the overall linear polarization in directions that are incompatible 
with the observed Q-U time series. 

Computing the model predictions relies on ray tracing, which 
requires knowledge of the path rays take in 3D curved space-time. These 
geodesic paths depend onthe unknown black-hole properties”: mass, 
spin and inclination. Nevertheless, the mass of Sgr A is well constrained 
through stellar dynamics”; M ~ 4 x 10° M,, where M, denotes solar 
mass. Furthermore, Fig. 1 (blue curve) illustrates that the data fit is not 
very sensitive to black-hole spin: a € [0, 1]. Thus, the only remaining 
unknown is the inclination angle. 

To estimate the inclination, we numerically bin 0, € [0, 1/2] and 
recover the 3D emission for every given (fixed) angle. For each angle, 


we recover a (locally) optimal 3D emission by minimizing a x? loss over 
the model parameters. Practically, for numerical stability, we avoid 
the extreme angles of face-on and edge-on by gridding 0, € [4°, 80°] 
(at 2° increments). Figure 1 plots the validation-x7: a robust likelihood 
approximation for 0,, which appears to favour low inclination angles 
(details in Methods). For each inclination, the recovery is run five times 
with a random initialization for the 3D structure. Therefore, the error 
bars are not a measure of posterior uncertainty; rather, they indicate 
the stability of the locally optimal solution. 

An overview of the tomographic data-fitting framework is illus- 
trated in Fig. 3. Mathematically, the 3D emission volume is estimated 
by minimizing a x” loss between the observed LP and the model pre- 
diction. The continuous 3D emission volume is represented by a 
coordinate-based, fully connected, neural network (‘neural represen- 
tation’) andis constrained toa domain witha radius of 6M < r < 20Mand 
close to the equatorial disk |z| < 4M (6M is the ISCO of a non-spinning 
black hole). The data fit used in this work relies on the reduced x defini- 
tions of ref. 29. This is not a strict definition of reduced y that includes 
a normalization by the degrees of freedom. Rather, it is normalized 
by the total number of data points, which is useful for comparing fit 
quality in our experiments where degrees of freedom remain fixed. 

The ill-posed inverse problem we solve does not have a unique 
solution. The recovered 3D structure depends on, among other factors, 
the assumed inclination angle. Furthermore, solving anon-convex opti- 
mization problem with stochastic gradient descent methods leads toa 
local (and not global) minimum. Thus, the recovered 3D structure also 
depends onthe random initialization of the network weights. Figure 4 
highlights the robustness of the recovered 3D structure across different 
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Fig. 3 | An overview of the orbital tomography framework based on light 
curve observations. (1) An orbital model propagates an initial (canonical) 
emission volume (e) in time. (2) Ray tracing: we compute GR ray paths according 
to the black-hole parameters and numerically integrate the 3D emissivities to 


synthesize image-plane frames. (3) Each frame is summed to produce a single 
light curve data point that, downstream, is compared to the observations. (4) A 
neural representation of the underlying 3D volume. Each componentis discussed 
extensively in Methods. 


inclinations and initial conditions. Although this is not an exploration 
of the posterior distribution, the different recoveries give a sense of the 
solution’s stability. Qualitatively, the details of each recovered struc- 
ture exhibit dependence on both the inclination angle and initializa- 
tion. Nevertheless, some key features are consistent across these two 
axes. Although the exact angular extent of the structures is not stable, 
the azimuthal and radial positions appear stable and consistent with the 
average structure. Moreover, the separation of the emission into two 
distinct structures appears consistent across the different recoveries. 

To analyse the ability to recover and detect different underlying 
3D morphologies, we simulated synthetic datasets mimicking ALMA 
observations for three underlying 3D structures: simple hotspot, flux 
tube, double source. Figure 5 highlights the recovery results obtained 
from these datasets at two (unknown) inclination angles. A compre- 
hensive analysis of the simulated datasets and reconstruction results 
is given in Supplementary Information Section 2. 

Recovering the 3D structure from light curve observations is 
highly ill-posed. Thus, the recovery relies on physical constraints and 
model choices that we impose through the gravitational and synchro- 
tron emission models. The robustness of the results depends on the 
validity of these model choices, detailed in Table 1 and discussed below. 

The key assumption for orbital tomography is that the emission 
is in orbit within an accretion disk near the equatorial plane and can 
be modelled as a simple transformation to a canonical (or initial) 3D 
emission. Note that small shifts from the equatorial plane are allowed 
by our model. This enables the formulation of an inverse problem for 
estimating the 3D emission from observations. We consider orbits 
characterized by a Keplerian angular velocity profile (neglecting radial 
or vertical velocity components), accounting for shearing due to differ- 
ential rotation (ignored by the previous analyses"”>”°) while neglecting 
the dynamics of cooling, heating, expansion and turbulence. Although 
this simplifying assumption does not hold in general, it is consistent 
with theoretical simulations’, which show consistent structures on 
short -1 orbit timescales. 

Furthermore, in modelling synchrotron emission, we assume a 
homogeneous vertical magnetic field that is externally fixed and is 
independent of the flare or accretion disk dynamics. The choice of a 
vertical magnetic field for the fiducial recovery is motivated by the 
notion that vertical magnetic fields could be powering Sgr A’ flares, 
apparentin GR magnetohydrodynamic simulations that produce mag- 
netic eruption events”. Moreover, from an observational standpoint, 
vertical magnetic fields are preferred by both the near-infrared analysis 
of GRAVITY"? and millimetre ALMA analysis of ref. 13. Nevertheless, 


Random initialization (seed) 
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Fig. 4 | A visualization of 3D recoveries across different inclinations and 
initial conditions. Although some of the details of the recovered structure 
depend on both axes, some key features remain consistent. The exact angular 
extent of the structures is not stable; nevertheless, the azimuthal and radial 
positions appear stable and consistent with the average structure highlighted 
in Fig. 1. Moreover, the separation into distinct emission regions of an elongated 
feature trailed by asmaller, dimmer, compact bright spot appears consistent 
across the different recoveries. 


the true spatial structure and dynamic properties of the magnetic fields 
around Sgr A are largely unknown. 

In Fig. 6, we analyse some of the systematic model choices detailed 
above by exploring the effects of (1) magnetic field configuration, (2) 
rotation direction and (3) sub-Keplerian orbits on the data fit and 3D 
reconstruction. It is important to note that we do not aim to exhaus- 
tively test all possible magnetic field and orbital velocity models, but 
instead highlight the sensitivity of our reconstruction to these model 
choices. The top-left panel of Fig. 6 compares the validation-y? for three 
magnetic field configurations: vertical, radial and toroidal, respec- 
tively, denoted by subscripts z, rand @. For aradial magnetic field, the 
best-fit recovery is not acompact bright emission region (Fig. 6, bottom 
left). Rather, it is a fainter, diffuse structure. Even so, according to the 
data fit and consistent with prior studies, vertical magnetic fields are 
preferred with a lower validation-x’ value. 

The centre and rightmost panels in Fig. 6 highlight how clockwise 
rotation (CW) and a Keplerian orbit are favourable to anticlockwise 
rotation (CCW) or sub-Keplerian orbits, consistent with the analyses of 
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Fig. 5 | 3D recoveries for three simulated structures observed at two 
(unknown) inclination angles, glo = 12°and phi = 64°. Using synthetically 
generated light curves as observations, the 3D reconstructions are able to recover 
different flare morphologies in the presence of background accretion noise (not 
visualized in this figure). Further analysis and details are given in Supplementary 
Information Section 2.3. 


GRAVITY” and ref. 13. To test the fit of orbit direction and sub-Keplerian 
fraction, we set the angular velocity profile to Q = +f,Q,, where the 
+ sign dictates the direction (CW/CCW) and f, the magnitude (f, =1 
results in a clockwise Keplerian orbit). The 3D recovery under f, = 0.9 
is shown in the bottom-right panel of Fig. 6, highlighting a broadly 
consistent recovery with the fiducial model assumptions. However, 
this consistency eventually breaks down at strong deviations fromthe 
Keplarian velocity assumption (also resulting in lower validation-x? 
values). The top view illustrates how a sub-Keplerian orbit impacts the 
recovered flare’s radial position. 

Another key assumption made by our work is that the millimetre 
emission region is optically thin. This is consistent with both EHT 
observations of Sgr A’ (ref. 17) and their theoretical interpretation”. 
Moreover, theoretical analysis has shown that a ‘flux tube’ flare would 
be optically thinner due to its higher temperature and lower density 
compared to the surrounding accretion flow. 

Finally, we assume that the recovered emission structure is within 
the accretion disk. Although our model does not account for an alterna- 
tive jet interpretation, the assumption of accretion flares is consistent 
with (1) theoretical simulations showing powerful equatorial current 
sheaths associated with flux eruptions forming within the accretion 
disk”, (2) observational evidence from EHT/very-long-baseline inter- 
ferometry analyses consistent with a compact source model without 
any detectable jet contribution”*°” and (3) the GRAVITY detection, 
indicating an astrometric centre aligning with the mass centre, which 
implies that the orbiting feature is in proximity to the equatorial plane. 
The alternative scenario demands a precise alignment between the 
direction of the jet and the observer's line of sight. 


Discussion 

We present a computational approach to image dynamic 3D struc- 
tures orbiting the most massive objects in the universe. Integrating 
polarimetric general relativistic ray tracing and neural radiance fields 
enables resolving a highly ill-posed tomography in the extremely 
curved space-time induced by black holes. Applying this approach 
to ALMA observations of Sgr A’ reveals a 3D structure of a flare, witha 


Table 1 | Summary of the key physical assumptions made in 
the modelling 


Emission model Synchrotron fixed vertical magnetic field; optically 


thin disk 


Dynamical model Keplerian; to=9:20 UT; clockwise orbit (no radial/ 


vertical velocity); velocity shear 


Gravitational model Kerr; mass=4.154x10°M,; non-spinning; 8, 


estimated from data 


3D model Neural representation; recovery domain: 6M <r<20M 


(FOV=200 arcsec); |z|< 4M 


We assume that the emission source is in orbit around a black hole within its accretion disk. 
The recovered 3D emission relies on constraining the flexibility of 3D neural fields with 
black-hole physics. Thus, the accuracy of the reconstruction depends on the validity of the 
model assumptions. Figure 6 explores the effects of some of the assumptions (magnetic field 
configurations, orbit direction and sub-Keplerian orbits) on both the data fit and recovered 
3D. We assume a non-spinning black hole because our analysis found that results are only 
weakly sensitive to black-hole spin (Fig. 1a). 


location broadly consistent with the qualitative analysis presented in 
ref. 13. This attempt at a 3D reconstruction ofa Sgr A flare suggests an 
azimuthally elongated bright structure at a distance of 11M trailed by 
a dimmer source at 13M. Although the recovered 3D is subject to and 
sometimes sensitive to the gravitational and emission models, under 
physically motivated choices, we find that the 3D reconstructions are 
stable and our approachis successful on simulated data. Moreover, our 
data-fit metrics provide constraints favouring low inclination angles 
and clockwise rotation of the orbital plane, supporting the analyses of 
ref. 13, EHT and GRAVITY”. 

Orbital polarimetric tomography shows great promise for 3D 
reconstructions of the dynamic environment around a black hole. 
Excitingly, extending the approach and analysis to spatially resolved 
observations (for example, EHT) and multifrequency data could ena- 
ble relaxing assumptions to further constrain the underlying physi- 
cal structures that govern the black hole and plasma dynamics (for 
example, black-hole spin, orbit dynamics, magnetic fields). To that 
end, future work will likely need to extend our model to non-optically 
thin media and non-azimuthal velocity patterns. Lastly, by adapting 
orbital polarimetric tomography to other rich sources of black-hole 
time series observations (for example, quasars and microquasars), 
this imaging technology could open the door to population statistics 
and improve our understanding of black holes and their accretion 
processes. 


Methods 
Inthe following section, we describe our methodology, which is evalu- 
ated on synthetic simulations and analysed in the Supplementary 
Information. 


Preprocessing 

We reduced the -100 min of ALMA light curves by time-averaging 
over ~1 min intervals, resulting in ~100 data points for each Stokes 
component. Following the procedures outlined in ref. 13, we subtract a 
constant LP component with magnitude and angle of Pjisx = 0.16 Jy and 
Eais = 37°, respectively, to account for the background accretion disk; 
we derotate the electric vector polarization angle by 32.2° to account 
for the estimated Faraday rotation’. We model the data as homosce- 
datic within a short and stable observation window witha polarimetric 
noise level estimated at 0g = 0, = 0.01 Jy (ref. 13). Although we do not 
fit the total intensity, we regularize the model to have a total intensity 
around 0.3 Jy with a standard deviation of 0.15 Jy (ref. 13). Following 
ref. 13, we set 9:20 UT as the initial time of the analysis and 3D recon- 
struction of the flare. Supplementary Fig. 3 shows an analysis of differ- 
ent initial times around 9:20 UT, which provides further motivation for 
the selection of this initial time. 
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Fig. 6 | The effects of different model choices: magnetic field, rotation 
direction and orbital velocity. Top: validation-y” under each model choice. 
Bottom: 3D reconstructions under various model assumptions. The red curve in 
all panels represents the fiducial model parameters: vertical magnetic field (B,), 
CW rotation and Keplerian orbit (fx = 1.0). The global minimum for each curve 

is highlighted by a horizontal dashed line in the respective colour. Left: three 
different magnetic field configurations: vertical, radial and toroidal (subscripts 
z, rand @, respectively). The recovered 3D under a radial magnetic field appears 
spread out rather than as a compact hotspot-like structure. That being said, 


consistent with the analysis of refs. 11,13, vertical magnetic fields, which do result 
inacompact hotspot-like structure, are favourable according to this metric, with 
lower validation-y? around 0, = 12°. Centre: a comparison of CW and CCW angular 
velocity models. Consistent with the analysis of ref. 13, a CW rotation is preferred 
across all inclination angles. Right: a Keplerian orbit has the lowest validation-x” 
fit across three different fractions of sub-Keplerian orbit: fx = 1.0, 0.9 and 0.8. The 
recovery under fx = 0.9 (bottom right) is broadly consistent with the Keplerian 
model, with a tendency towards smaller radii (illustrated by the top-view panel). 


Forward model 

Inthis section, we formulate the forward model, which takes a canoni- 
cal 3D emission around a black hole as input and synthesizes light 
curves as output. Figure 3 provides a high-level overview of the forward 
model divided into four key building blocks, which we describe in the 
sections below. 


Orbit dynamics. The key assumption for orbital tomography is that 
the four-dimensional emission e(t, x), where t denotes time and x 
denotes 3D spatial coordinates,is in orbit around the black hole and 
can be modelled as a simple transformation of a canonical (or initial) 
3D emission, €)(x): 


e(t, x) = €o(T:X). (1) 


The transformation T, propagates the initial 3D structure in time and 
connects dynamic observations, such as light curves, to the canoni- 
cal 3D structure. This in turn enables formulating an inverse problem 
of estimating e,(x) from time-variable observations. Although the 
assumption of a coordinate transformation does not hold in general, 
it is well suited for compact, bright structures over short time scales, 
during which complex dynamics are negligible. 

In our work, we consider a Keplerian orbit model with an angular 
velocity 


VM 


20) = —" _, 
I2 + ayM 


(2) 


whereris the distance from the black-hole centre and Mis the black-hole 
mass. Note that for a = O, equation (2) coincides with the Newtonian 
expression for angular velocity. A purely azimuthal orbit is suitable 


outside the ISCO, where radial velocities play a smaller role. Thus, we 
formulate the coordinate transformation as a shearing operation: 


T,= So: (3) 


where S, is a rotation matrix at an angle 


(6r) = (C= bo) 200). (4) 


The angular velocity dependence onr (equation (2)) causes shearing 
due to the faster motion of inner radii. 


Image formation. In this section, we describe how e, relates to light 
curve observations through an image-formation model. Each image 
pixel collects radiation along a geodesic curve: /(0, a, 8) terminating at 
the image coordinates (a, f). The ray path Tis determined by a handful 
of black-hole parameters: © = (M,a, 6,). Omitting the explicit depend- 
ency onimage coordinates (for brevity), we model image pixels through 
the polarized radiative transfer” ** ofan optically thin disk (attenuation 
can be neglected for Sgr A’ 230 GHz observations”): 


pit) 
Patt) 
p) = = 
Pu) 
pv(o) 


gx elt + T, X)R(x)J(x)dx. (5) 


xel(0) 


Equation (5) describes how pixel values are computed through an 
integration along geodesic curves / computed by solving a set of dif- 
ferential equations” (Supplementary Information Section 3). The 
integrand comprises four elements: g, e, R and J. Here e(t + Tą, X) isthe 
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unknown scalar emissivity that depends on microphysical properties 
(for example, local electron density and temperature) and txis the time 
delay that accounts for photon travel time (often referred to as slow 
light). We model the polarized synchrotron radiation as this scalar 
emissivity function multiplied by a Stokes-vector J proportional to” 


Ji «8 (\B| sin bp)” (6) 
Ja «OJ; (7) 
Jy =0 (8) 


In this work, we consider only linear polarization, thus setting J, = 0. 
Moreover, the spectral index (reflecting the change in the local emis- 
sion with frequency) is approximated as a, ~ 1 (ref. 36). Note that the 
local emission frame is defined to align with Stokes Q; therefore, Jy = 0. 
Thescaling factor q;€ [0, 1] is the (volumetric) fraction of linear polari- 
zation, and @, is the angle between the local magnetic field B and pho- 
ton momentum k, given by 


k(x) x B(x) 


soba: OOE 9 
Ik@O|1B00| 0) 


sin h(x) = 


The two remaining quantities to define are R and g. The matrix R 
rotates the LP, (Jo, Ju), from the emission frame to the image coordi- 
nates through parallel transport” (Supplementary Information Sec- 
tion 3.4). The scalar field g(x) is a GR red-shift factor, which decreases 
the emission when the material is deep in the gravitational field or 
moving away from the observer. More generally, g(x) depends on the 
local direction of motion, u, relative to the photon momentum k: 


&(X) = ( u(x), k09) (10) 


Note that u and kare 4-vectors, more explicitly defined in Supplemen- 
tary Information Section 3. 


Light curves. For a given 3D emission, ray-tracing equation (5) enables 
computing a single pixel value over time. We compute light curves by 
numerically sampling a large field-of-view (FOV) and summing over 
image-plane coordinates: 


10) 
IQ) 
IQ = = 2, P(t, a, p) (11) 
Iy(t) 2 
I(t) 


Neural representation. We formulate a tomographic recovery relying 
onaneural representation” of the unknown 3D volume: e€,(x). Thus, 
instead of a traditional voxel discretization, the volume is represented 
by the weights, w, of a multilayer perceptron (MLP), that are adjusted 
to fit the observations. 

The implicit regularization of the MLP architecture enables tack- 
ling highly ill-posed inverse problems””*’. The MLP takes continuously 
valued coordinates x as input and outputs the corresponding scalar 
emission at that coordinate 


€o (X) = MLP y(V(x)), (12) 


where y(x) is a positional encoding of the input coordinates. 

Studies have shown” that encoding the coordinates instead 
of directly taking them as inputs can capture continuous fields 
better (converging in the width limit to a stationary interpolation 
kernel”). Thus, our work relies on a positional encoding that projects 


each coordinate onto aset of sinusoids with exponentially increasing 
frequencies: 

y(x) = [sin(x), cos(x), ..., sin (2171x) , cos (2!-x)] (13) 
The positional encoding controls the underlying interpolation kernel 
used by the MLP, where the parameter L determines the bandwidth of 
the interpolation kernel”. 

In our work, we use a small MLP with four fully connected layers, 
where each layer is 128 units wide and uses ReLU activations. We usea 
maximum positional encoding degree of L = 3. The low degree of L is 
suitable for volumetric emission fields, which are naturally smooth”. 

Once the neural network weights, w, are adjusted to fit the data, 
the network can be sampled at any 3D point, x, to produce the emis- 
sion value at that point. This allows us to sample the network at regular 
grid points to extract a 3D volume representation of the recovered 
emission. Moreover, we can sample the network along straight-ray 
paths and ray trace the recovered emission as it would be seen by a 
perspective (pinhole) camera in flat space (used for the 3D visualiza- 
tions throughout the Article). 


Solving the inverse problem 

In this section, we formulate an optimization approach that enables 
jointly estimating the 3D emission and inclination, which are the param- 
eters of the forward model. Supplementary Fig. 1 shows a high-level 
illustration of the data-fitting procedure introduced in the following 
section. 


Tomographic reconstruction. To estimate the 3D emission from light 
curve observations, we formulate a minimization problem. We estimate 
w, which parameterizes e,(x), by minimizing a x’ data fit for each Stokes 
component, evaluated for a fixed set of black-hole parameters, ©: 


X? (WIO) = xX; (WIO) + xz (WI) + x% (w10). (14) 


Here, we restrict the discussion to the total intensity and LP compo- 
nents: /, Q, U. Each x” is calculated as a sum over discrete temporal 
data points 


2 
EAL as 


1 
X¢(wlO) = Nos a J; 


obs 


where N,,, is the total number of data points; the subscript s = {/, Q, U} 
represents the stokes components; and y,, /,and o, are the polarimetric 
observations, model and noise standard deviation, respectively. Note 
that /,(¢, w|O) is simply the light curve given by equation (11), sampled 
at discrete time t, where w|O highlight its dependency/conditioning 
onthe network/black-hole parameters. 

Equation (14) depends on unknown black-hole parameters; 
nevertheless, the mass of Sgr A’ can be constrained through stellar 
dynamics”®*°; M = 4 x 10° M,, where M, denotes solar masses. Further- 
more, because the data fit is insensitive to black-hole spin, the only 
estimated parameter is the inclination angle. To estimate the inclina- 
tion, we numerically bin 6, € [0, 1/2] and recover the 3D emission by 
minimizing equation (14): 


w* (0o) = arg min x? (w|9,) . (16) 


By interpreting equation (16) as a function of 0,, we approximate the 
marginal log-likelihood as 


L(Ooly) x X?(Oo|w*). (17) 


Equation (17) isa zero-order expansion about the maximum likelihood 
estimator: w*. 
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Model selection using validation-y*. Although equation (17) 
tells us how well each model (inclination) fits the data, it is susceptible 
to overfitting. To mitigate overfitting, we define a more robust metric 
called validation-x” (Supplementary Fig. 5). The inclination angle is 
then estimated through the following procedure: 


(1) During optimization, ray positions are fixed to the centre of 
each image pixel. In our recoveries, we use an evenly sampled 
64 x 64 grid for a FOV of 200 parcsec. 

We compute y’ for perturbed pixel positions (off-centre) within 
a small pixel area. In our recoveries, we used a pixel area of 
3.125 x 3.125 arcsec’. 

We average x° of 10 randomly sampled (uniform) ray positions 
to compute validation-x’ curves. 

6* is estimated as the global minimum of the validation-x’. 


(2) 


(3) 
(4) 


Through simulations, we highlight how this procedure is amore 
robust selection criterion for models that are not overfitting the fixed 
ray positions (Supplementary Fig. 5). 


Optimization procedure. The neural network was implemented in 
JAX“. Both the synthetic experiments (Supplementary Information 
Section 2) and the ALMA recovery were optimized using an ADAM 
optimizer” with a polynomial learning rate transitioning from 
1x10*>1x10°° over 50,000 iterations. Run times were -1 h on two 
NVIDIA Titan RTX GPUs. Network weights were randomly initialized 
(Gaussian distributed) with several initial seeds. 


Data availability 

This paper makes use of the ALMA dataset ADS/JAO.ALMA# 
2016.1.01404.V, available through the ALMA data portal. Fully cali- 
brated data and other materials are available from the corresponding 
author upon reasonable request. 


Code availability 

The software packages used to analyse the data are available at the 
following sites: kgeo (https://github.com/achael/kgeo) and bhnerf 
(https://github.com/aviadlevis/bhnerf). 
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